clear all
set mem 100m
local path0="C:\Dropbox\GeneticsProject (1)\REStat\round_accepted_replication_files\"
local path1="C:\Dropbox\GeneticsProject (1)\REStat\round_accepted_replication_files\outregs"

cd "`path0'Workfiles\"
set more off
use completed_data, clear

xi: reg  pwt_ln_rgdpwok idv, robust
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		replace bdec(3) rdec(3) aster nocons nolabel ctitle("OLS") 
		

xi: ivreg2  pwt_ln_rgdpwok (idv= distM_UK), first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + distance") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
		
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*							direct gene effect
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* exclude small countries (here Singapore)
xi: ivreg2  pwt_ln_rgdpwok (idv= HTTLPR) if countrycode~="TTO", first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + httlpr_s") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
	
xi: ivreg2  pwt_ln_rgdpwok (idv= HTTLPR distM_UK) if countrycode~="TTO", first gmm2s robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + httlpr_s") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
	
	local J_httlpr_s=e(j)
	
*~~~~~~~~~~~~~~~~~~~~~~~~~~		
* exclude outlier (Nigeria)
* results are similar when huber robust regression is used in the first stage
xi: ivreg2  pwt_ln_rgdpwok (idv= A118G ) if countrycode~="NGA", first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + allele_g") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
	
xi: ivreg2  pwt_ln_rgdpwok (idv= A118G distM_UK) if countrycode~="NGA", first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + allele_g") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')

	local J_allele_g=e(j)
				
	
*~~~~~~~~~~~~~~~~~~~~~~~~~~		
xi: ivreg2  pwt_ln_rgdpwok (idv= pathogen_ms_9_items), first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + pathogen_ms_9_items") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
		
xi: ivreg2  pwt_ln_rgdpwok (idv= pathogen_ms_9_items distM_UK) if countrycode~="SGP", first robust
	matrix A=e(first)
	local R2first=A[2,1]
	local Ffirst =A[3,1]	
	outreg2 idv using "`path1'\RR_table001.txt" , ///
		append bdec(3) rdec(3) aster nocons nolabel ctitle("IV + pathogen_ms_9_items") ///
		addstat("1st stage F-stat", `Ffirst',"Partial R2",`R2first')
		
	local J_pathogen_contemporary=e(j)

	
	
*=============================================================================
*			This part generates critical values for the J-test
*			using boostrap with recentering as in Hall and Horowitz
*=============================================================================

use completed_data, clear
set seed 1234567 

* exclude outlier: Trinidad & Tobago 
ivreg2  pwt_ln_rgdpwok (idv= HTTLPR distM_UK) if countrycode~="TTO", first gmm2s robust
local NN=e(N)
keep if e(sample)
gen z1=HTTLPR
gen z2=distM_UK
local J_httlpr_s=e(j)

keep pwt_ln_rgdpwok idv z1 z2
save temp_book0, replace

ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust

matrix eW=e(W)
predict res0, resid

gen m0=z1*res0
sum m0
local m0s=r(mean)

gen m1=z2*res0
sum m1
local m1s=r(mean)

gen m2=res0
sum m2
local m2s=r(mean)

matrix A0=(`m0s' \ `m1s' \ `m2s')

matrix JJ=A0'*eW*A0*`NN'
local J0=JJ[1,1]
di `J0'

**** start the bootstrap cycle

tempname sim
postfile `sim' JX JY using boot_results_httlpr_s, replace
quietly {
	forvalues i = 1/1000 {
		use temp_book0, clear
		bsample
		ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust
		
		matrix eW=e(W)
		matrix eS=e(S)
		predict res0, resid

		gen m0=z1*res0
		sum m0
		local m0s=r(mean)

		gen m1=z2*res0
		sum m1
		local m1s=r(mean)

		gen m2=res0
		sum m2
		local m2s=r(mean)

		matrix A1=(`m0s' \ `m1s' \ `m2s')


		matrix JJ=(A1-A0)'*eW*(A1-A0)*`NN'
		
		matrix Sc=eS-A0*A0'-A1*A0'-A0*A1'
		matrix eWc=inv(Sc)
		matrix JJc=(A1-A0)'*eWc*(A1-A0)*`NN'
		
		local JX=JJ[1,1]
		local JY=JJc[1,1]
		
		if mod(`i',100)==0 {
			noisily di "`i'"
			}
			
		post `sim' (`JX')  (`JY')
	}
}
postclose `sim'

*=========================================================================


use completed_data, clear
set seed 1234567 


* exclude outlier: Nigeria
ivreg2  pwt_ln_rgdpwok (idv= A118G distM_UK) if countrycode~="NGA", first gmm2s robust
local NN=e(N)
keep if e(sample)
gen z1=A118G
gen z2=distM_UK
local J_allele_g=e(j)

keep pwt_ln_rgdpwok idv z1 z2
save temp_book0, replace

ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust

matrix eW=e(W)
predict res0, resid

gen m0=z1*res0
sum m0
local m0s=r(mean)

gen m1=z2*res0
sum m1
local m1s=r(mean)

gen m2=res0
sum m2
local m2s=r(mean)

matrix A0=(`m0s' \ `m1s' \ `m2s')

matrix JJ=A0'*eW*A0*`NN'
local J0=JJ[1,1]
di `J0'

**** start the bootstrap cycle

tempname sim
postfile `sim' JX using boot_results_allele_g, replace
quietly {
	forvalues i = 1/1000 {
		use temp_book0, clear
		bsample
		ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust
		
		matrix eW=e(W)
		predict res0, resid

		gen m0=z1*res0
		sum m0
		local m0s=r(mean)

		gen m1=z2*res0
		sum m1
		local m1s=r(mean)

		gen m2=res0
		sum m2
		local m2s=r(mean)

		matrix A1=(`m0s' \ `m1s' \ `m2s')


		matrix JJ=(A1-A0)'*eW*(A1-A0)*`NN'
		local JX=JJ[1,1]		
		if mod(`i',100)==0 {
			noisily di "`i'"
			}
			
		post `sim' (`JX')
	}
}
postclose `sim'



*=========================================================================
*					historical pathogen prevalence
*=========================================================================

use completed_data, clear
set seed 1234567 

ivreg2  pwt_ln_rgdpwok (idv= pathogen_ms_9_items distM_UK), first gmm2s robust
local NN=e(N)
keep if e(sample)
gen z1=pathogen_ms_9_items
gen z2=distM_UK

keep pwt_ln_rgdpwok idv z1 z2
save temp_book0, replace

ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust

matrix eW=e(W)
predict res0, resid

gen m0=z1*res0
sum m0
local m0s=r(mean)

gen m1=z2*res0
sum m1
local m1s=r(mean)

gen m2=res0
sum m2
local m2s=r(mean)

matrix A0=(`m0s' \ `m1s' \ `m2s')

matrix JJ=A0'*eW*A0*`NN'
local J0=JJ[1,1]
di `J0'

**** start the bootstrap cycle

tempname sim
postfile `sim' JX using boot_results_path_cont, replace
quietly {
	forvalues i = 1/1000 {
		use temp_book0, clear
		bsample
		ivreg2  pwt_ln_rgdpwok (idv= z1 z2), first gmm2s robust
		
		matrix eW=e(W)
		predict res0, resid

		gen m0=z1*res0
		sum m0
		local m0s=r(mean)

		gen m1=z2*res0
		sum m1
		local m1s=r(mean)

		gen m2=res0
		sum m2
		local m2s=r(mean)

		matrix A1=(`m0s' \ `m1s' \ `m2s')


		matrix JJ=(A1-A0)'*eW*(A1-A0)*`NN'
		local JX=JJ[1,1]		
		if mod(`i',100)==0 {
			noisily di "`i'"
			}
			
		post `sim' (`JX')
	}
}
postclose `sim'



**=================================================================
*					report p-values
**=================================================================
noisily di "httlpr_s"
use boot_results_httlpr_s, clear
gen t0=(JX>`J_httlpr_s')
sum t0

noisily di "allele_g"
use boot_results_allele_g, clear
gen t0=(JX>`J_allele_g')
sum t0

noisily di "Historical pathogens: extended"
use boot_results_path_cont, clear
gen t0=(JX>`J_pathogen_contemporary')
sum t0	
